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Abstract 

Assuming a non-paraxial propagation operator, we study the propagation of an electromagnetic field with 
an arbitrary initial condition in a quadratic GRIN medium. We show that at certain specific periodic distances, 
the propagated field is given by the fractional Fourier transform of a superposition of the initial field and of a 
reflected version of it. We also prove that for particular wavelengths, there is a revival and a splitting of the 
initial field. We apply this results, first to an initial field given by a Bessel function and show that it splits into 
two generalized Bessel functions, and second, to an Airy function. In both cases our results are compared with 
the exact numerical ones. 


1 Introduction 


Graded index (GRIN) media are mainly used in image formation applications [^. On the other hand, it has been 
established that a GRIN medium can form self-images of periodic fields 2][3 and can support invariant propagation 
modes, either in the paraxial and the non-paraxial domains |^. In a different context, light propagation in a 
quadratic GRIN medium can be employed as a form of optical emulation of quantum phenomena. An example is 
the mimicking of quantum mechanical invariants by the propagation of light through the interface of two coupled 
GRIN devices . Because the Schrodinger equation and the paraxial wave equation in classical optics are formally 
equivalent, cross applications between quantum mechanics and classical optics are common. One can extend the 
application in order to consider not only the paraxial regime but also the non-paraxial one, i.e. the complete 
Helmholtz equation. For instance, supersymmetric methods, common to quantum mechanics, have been proposed 
in classical optics [^|^. 

An interesting conceptual and mathematical result is that the paraxial propagated field in a quadratic GRIN 
medium can be expressed as the Fractional Fourier transform (FrFT) of the incoming beam [^. Indeed, the FrFT 
operator has been considered for discussing additional similitudes between classical optics and quantum formalism. 
For instance, Agarwal and Simon have shown that Fresnel diffraction leads to the FrFT by noting that, when 
constructing the (quantum) harmonic oscillator evolution operator, it contains a term proportional to paraxial free 
propagation. In another report, by Fan and Ghen 10 , it is shown that the quantum-mechanical position-momentum 


mutual transformation operator is the core element for constructing the integration kernel of FrFT. 

In the context of GRIN media, with quadratic refractive index dependence in the radial coordinate, knowledge 
from harmonic oscillator-type Hamiltonians can be used for the solution beyond the paraxial regime. As a significant 
result in this context, we established previously the long period revival and splitting of a Gaussian beam, transmitted 


by a quadratic GRIN medium II , 12 . Such effects are theoretically predicted expressing the Taylor series for the 


propagation operator up to the second order, i. e. including an additional term beyond the paraxial approximation. 

In the present paper, we describe a significant generalization of such long period effects, assuming again a second 
order form of the propagator. We establish propagation distances that exhibit the revival and the splitting of an 
arbitrary input field. We also establish propagation distances for which any input field E{x) splits into two fields 
given by the FrFT of E{x) + E{—x). We apply this result to an initial field given by a Bessel function, and as a 
solution, we obtain the superposition of two so-called Generalized Bessel functions IS -Tsl. 
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2 Helmholtz equation for GRIN media 

The Helmholtz equation in two dimensions for a GRIN medium is 


■ 52 92 




EiX,Z)=0, 


( 1 ) 


where k is the wave number and n{X) is the variable refraction index. For a quadratic medium, the refraction 
index can be written as 

n^{X) = nl{l-g^X^), ( 2 ) 

where g is the gradient index in the X direction. So, for a quadratic dependence in the index of refraction, the 
Helmholtz equation is expressed as 


^ ^ ^) = 


(3) 


where we have introduced the dimensionless variables x = ^/gX and z = and where we have defined g = nogn 
and k = tlqk/ y/g. 

Introducing the operator p = —ij-, we can cast the last expression as 


19 


d F/(x, z) / 2 ^2 2 \ 7 - 1 / \ 

=-{k^-p^- X^) E{x, z), 


(4) 


whose formal solution is 




P^ + x'^ 

—iz\/k‘^ — fZ — x"^ 

Eix, 0 ) = exp 



E{x,G), 


(5) 


where E{x, 0) is the boundary condition for z = 0. 

We define the lowering ladder operator, a = (1/2)^/^ (a; + fp), its adjoint, the raising ladder operator, a'^ = 
(1/2)^/^ {x — ip) and the number operator n = a^a, and we write Eq. ([^ as 


E(x, z) = exp 


—izk 


2h + 1 


E{x,Q), 


( 6 ) 


2.1 The paraxial approximation. 

As a background for our main result in next section, we present here an alternate derivation of the paraxial 


propagation in GRIN media, in terms of the FrFT 10 21 . The paraxial approximation is obtained when the square 


root in the exponential of Eq. ([^ is expanded to first order, that is 


E{x, z) = exp 


—izk 1 — 


2/c2 


exp E(x, 0 ). 


(7) 


On the other hand, it has been established that the fractional Fourier transform of a well behaved function f{x) 
can be obtained in terms of the number operator h, as S'q {fix)} = exp (iah) {/(x)}, where alpha is the transform 
order. Considering this result, Q can be rewritten as 


E{x,z) = exp 






( 8 ) 


Thus, the paraxial propagation to a distance z is proportional to the fractional Fourier transform of order ^ of the 
initial condition E(a;, 0 ) (W 21 . 
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2.2 Beyond the paraxial approximation 

We now allow ourselves to go one step further than the paraxial approximation. In Eq. ([^, we again expand the 
square root in Taylor series, but we hold terms to second order instead, to obtain 


/ .71 ^ 

(.12 A 


■ 1 - 2 



exp 



E{x,0), 


where, for simplicity, we have defined 


7 4 1 7 2 1 


We develop the initial condition in terms of the Gauss-Hermite functions 




1/4 


where (x) are the Hermite polynomials, to obtain 


exp ( --x"^ ) (x), 


3=0 


E{x, z) = exp exp (^ij^znj ^ Cj exp 

Now, for z = Ink^ with I any non-negative integer, we have 
E{x, z = Ink^) = exp (— ilqiTr) exp n) 




tpj{x). 


^ C2,‘P23{x) +i^Yl C2j + l‘P2j + l{x) 
3=0 3=0 


(9) 


( 10 ) 


( 11 ) 


( 12 ) 


(13) 


Next, considering the identities 

OO OO . 

2j^C2j(p2j(x) =J^Cjipj(x) + (-iyJ^Cj(pj(x) = E(x,0) + E(-x,0) 

and 


3=0 


3=0 


3=0 


we obtain 


2j^C2j+i(p2j+i(x) = J^Cj(pj(x) - (-iyJ2cj(pj(x) = E(x,0) - E(-x,0); 

3=0 3=0 3=0 


E{x, z = lirk^) = exp (—ilqiTr) [—^— exp {il^ 2 'xn) E(x, 0) H- - — exp (4/7271 n) E{—x, 0)]. 


But, as we already said, {f{x)} = exp (iah) {f{x)} [^, thus 

1 + i^ 


E{x, z = Inky = exp (—4/7171) 




72 7r 


l-i'- 

E{x,Q) + -^—E{-x,0) 


(14) 


(15) 


(16) 


(17) 


Hence, at these periodic distances the field is the fractional Fourier transform of a superposition of the initial field 
and its specular image. It is clear that if the initial condition is symmetric, E{—x,Q) = E{x,0), then at those 
periodic distances, we will have just the fractional Fourier transform of it. In particular, when / is congruent with 0 
modulo 4, the field is the fractional Fourier transform of the initial condition times a phase factor. If / is congruent 
with 1 or with 3 modulo 4, we get the fractional Fourier transform of a superposition of the initial field and its 
specular image. In the case of / congruent with 2 modulo 4, we obtain a phase factor times the fractional Fourier 
transform of the specular image of the initial condition; of course, in this last case, if the initial condition is sym¬ 
metric we will have just the fractional Fourier transform of the initial condition. 


An interesting case is obtained when the dimensionless wave vector k is chosen such that lTr"f 2 = 2m7r, where 
now m is another positive integer. As ^ 2 Tvm, is the identity operator, we get 


E{x, z = Inky = exp {—il^c'x) 


-E{x,0) 


1-4* 


E{-x,0) 


(18) 
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where = {2m/l — 1/2)^^^ and 7 c = 3/8 + m (4m — 31) /P. Thus, for those periodic distances and values of k, 
we can obtain the revival of the initial condition (when I is congruent with 0 mod 4), a superposition of the initial 
condition and its specular image (when I is congruent with 1 or with 3 mod 4) and the specular image of the 
initial condition (when I is congruent with 2 mod 4). A similar result occurs when lTrj 2 = (2m. + 1) tt, where again 
TO is a positive integer, for which the fractional Fourier transform becomes the parity operator; but in this case 
E{x,0), in Eq. ( [T^ is replaced by E{—x,0), and vice versa. The easiest situation is when we pick I = 1, and then, 

z = TT (2m — I) , k = J2m — ^ and 71 = to (4to — 3) + |. Below, we will study the Bessel functions and the 
Airy function as initial conditions, and this particular case will be exemplified. 


3 A Bessel function as initial condition 


In the particular case of a Bessel function as initial condition, E(x, 0) = J,^(bx + a), where v is non-negative integer, 
we know from the Appendix A its fractional Fourier transform, Eq (35), thus 


E (x,z = = exp i 


- ( 71 + ^ 72 ^ lTT+^fx^ + ^] tan (It:j 2 ) 


-y/sec (Ittj 2 ) X 


l + i^ 


X b sec (Ittj2) + ^ (It^I'2) ; —* 






—X b sec (111^2) + ^ (^ 7 J‘ 72 ) ; i 


(19) 


where is the second order generalized Bessel function, defined in Eq. (32) of Appendix A. 

In Eigure[T] we show the field intensity at the first splitting distance z = irk^, when the initial condition is a Bessel 


i^r 


i^r 




Figure 1: The square of the propagated electric field at the first critical distance when the initial condition is 
Jy (x + 3). The quadratic GRIN media parameters are no = 1-5 and g = 10mm“^. The beam has k = 1099.7. The 


black continuous line is the graphic of expression (19) and the red dotted line is the exact numerical solution. 


function Ji,(x 3). The parameters of the quadratic GRIN medium are no = 1.5 and g = 10mm and we have 


taken k = 1099.7. The black continuous line is the fractional Fourier transform given in Eq. (19) and the red dotted 
line is the exact numerical solution. 


For the special case indicated in Eq. (18), we have 


E(x, z = Ink^) = exp (—ilycTr) 


1 -|- 2^ 1 — Z 

—-— Jv(bx -I- a) -I--— Jv(-bx + a) 


( 20 ) 


In Eig. we plot (20) (black continuous line) and the exact numerical solution (red dotted line), when 6=1 and 
a = 3. The GRIN media parameters are the same as in Eig. and we took I = 1 and to = 8 x 10®, so kc = 1264.91. 
The result of the propagation is just a linear combination of the initial condition with its specular image. 
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Figure 2: The square of the propagated electric field at the first critical distance when the initial condition is 
Ju {x + 3). The quadratic GRIN media has no = 1.5 and g = 10mm“^. We have taken 1 = 1, m = 8 x 10® and 
then kc = 1264.91. The black continuous line is the graphic of expression (20) and the red dotted line is the exact 
numerical solution. 


4 An Airy function as initial condition 


We take now as initial condition the Airy function Ai{bx + a). Considering the fractional Fourier transform of this 
initial condition, Eq. (45) in Appendix B, the field propagated to the periodic distances z = Ink^ is given by 


E {x,z = hrk^) = ^/secJl^^^ exp < i 


11 6 ® 

“2^^72 + 2 (^^ 72 ) + tan® (iTrya) 


1 + i' 


■ exp 


6® 


+ - 


1-i' 


■ exp 


—tan ( 17172 ) sec (/7r72) Ai 

■ 5 ® 1 

7—a; tan ( 17172 ) sec (l7r72) Ai 


^4 

bx sec (l7r72) + a —— tan^ 

^4 

—bx sec (l7r72) + a —— tan^ 


( 21 ) 


In Figurej^ we show the field intensity at the first splitting distance 2 ; = Trfc®, when the initial condition is an Airy 



Figure 3: The square of the electric field at the critical distance when the initial condition is Ai (a:). The quadratic 
GRIN media parameters are uq = 1.5 and g = 10mm“^, and k = 1000.5. The black continuous line is the graphic 


of expression (21) and the red dotted line is the exact numerical solution. 


function Ai(a;). The parameters of the quadratic GRIN medium are no = 1.5 and g = 10mm and we have taken 


k = 1000.5. The black continuous line is the fractional Fourier transform given in Eq. (21) and the red dotted line 
is the exact numerical solution. 


For the special case indicated in Eq. (18), we have 
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Figure 4: The square of the propagated electric field at the critical distance when the initial condition is Ai(a:). 
The quadratic GRIN media has ng = 1.5 and g = 10mm“^. We have taken I — 1, m = Ax 10® and then k = 894.43. 
The black continuous line is the graphic of expression (22) and the red dotted line is the exact numerical solution. 


E{x, z = Inkl) 


exp 


1 + i^ 


Ai(a; + 1) 


l-i' 


Ai(—a; + 1) 


( 22 ) 


In Fig. 1^ we plot ( [^ (black continuous line) and the exact numerical solution (red dotted line), 
parameters are the same than in Fig. and we took I = 1 and m = 4 x 10®, so k = 894.43. 
propagation is just a linear combination of the initial condition with its specular image. 


The GRIN media 
The result of the 


5 Conclusions 


We have shown that the propagation of an initial held distribution in a quadratic GRIN media, beyond the paraxial 
wave approximation, produces the revival and the splitting of the input held, at specihc propagation distances. It 
is also proved that the held at another propagation distances is given by the fractional Fourier transform of the 
superposition of the initial held with a rehected version of it. We have applied the results to an initial Bessel 
function and found that the propagated held is given by a superposition of Generalized Bessel functions 13 14 


Also, the example when the initial held is an Airy function is examined. In both concrete cases, our predictions are 
checked against the exact numerical solution and good agreement has been established. 


A The fractional Fourier transform of a Bessel function of the first 
kind of integer order 


It is known 21 that the fractional Fourier transform Sa{f} of an arbitrary function / is given by, IJal/} = 
where h is the number operator of the quantum harmonic oscillator. 

To hnd the fractional Fourier transform of an integer order Bessel function of the hrst kind, we use the Jacobi-Angers 


integral representation 24 25 


Then we have. 


TT 

Jn{bx + a) = — / dTexp{i [nr — {bx -I- a) sin(r)]}. 

2Tr J 

— TT 

■ - 1 r 

Sa{Jn{bx + a)} = / dr exp{j [nr — (6x J-a) sin(T)]} 

1 r 

= — / dr exp (mr) exp {iah) exp [—i {bx + a) sin(r)]. 

27r J-TT 


(23) 


(24) 


We now write the number operator as?T.= ^-|-^— i, and so exp (ian) = exp (—if) exp [if + x^)]. But we 
know that 

exp [iC {f + x'^)] = exp [if{C)x‘^] exp [-ip(C) {xp + px)] exp [if (Of] , (25) 
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where /(C) = i tan(2C) and g(() = | In [cos(2C)]. Thus, 


Sa {Jn(bx + a)} = 


exp (-if) exp [if (f) x"] 


27r 


dr exp (inr) exp —ig (xp + px) 


exp 


exp [-* ( 6 x + a) sin(T)]. 


But, 


exp 


if ^ p^ exp [—i {bx + a) sin(r)] = exp if ^ sin^ (r) exp [—i (bx + a) sin(r)], 


then, we have 


ga {Jfiibx + a)} ^ 2 ) [^/ ( 2 ) ^ ] y cJr exp (inr) exp i/5^ sin^ (r) 


exp 




Also 


exp 


ig (I) (xp + px) exp [-i {bx + a) sin(r)] = exp -g exp |-i (^exp -2g 


bx + a] sin 




then the fractional Fourier transform of the Bessel functions of the first kind can be written as 


( .Oi\ 


r.^/ax 2 I 



f OL\ 

exp(-*-^ 

exp 

y(2)*l 

exp 

-g 

^ 2 ). 


X—y dr exp (inr) exp i/^ —^ 6 ^ sin^ (r) exp|—i^exp —2g(^—^ &a; + sin (r) | . 

Writing sin^ (r) = 5 [1 — cos(2r)] and changing the integration variable from r to —r in the last formula, we 
to 


^a{Jn{bx + a)} = exp “d 


exp 


.a 




+i.t \ ^ ] X 


I, 


dr exp {—inr) exp 


/ (^) ^^cos( 2 r) exp|i(^exp -2g fox + a) sin (r)|. 


Introducing now the Generalized Bessel Functions, defined as 13 


4 ™^ {x, y,c)= c'-Jn-mi (x) Ji {y ), 

I —— 00 

and its integral representation 

1 

(x, J/; e*®) = — / dip exp [ix sin p + iy sin {mp + 9) — inp] 
J — 00 


we can cast (311 as 


T?a {Jn{bx + a)} = exp 




7(2) 

U ^ 


e-29(f) 


/(i). 


bx + a,—f{-];-i 


Substituting the functions /(/) = | tan(2C) and g{f) = pn [cos(2C)], we finally obtain 


T?a {Jn{bx + o)} = -v/secQ- exp - 


2 b' 

—a + 1 a; I tana 


[bx sec a + o, — tan a; —i 


(26) 

(27) 

(28) 

, (29) 

(30) 
arrive 

(31) 

(32) 

(33) 

(34) 

(35) 
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B The fractional Fourier transform of an Airy function 


We use again the known fact 21 that the fractional Fourier transform Sa{f} of an arbitrary function / is given 
by, Sa{f} = where h is the number operator of the quantum harmonic oscillator. 

To find the fractional Fourier transform of an Airy function, we use the integral representation 24 25 

OO 

1 /* / \ 

Ai{bx -\- a) = — / dr exp [ir {bx + a)] exp ( z — j . 


(36) 


Then we have. 


Uc {Ai( 6 a; + a)} = [ dr exp {irbx) exp (ira) exp [ 


27r 


1 

27r 


dr exp ( *^ ) 6 xp (ira) e*“" exp (irbx). 


(37) 


We now write the number operator ash=^ + ^— i, and so exp (ian) = exp (—if) exp [if )]. But we 

know that 

exp [iC {f + a;^)] = exp [i/(C)a;^] exp [-ig{C) {xp + px)] exp [i/(C)p^] , (38) 

where /(C) = 5 tan(2C) and g(C) = 5 In [cos(2C)]. Thus, 

exp(-*f)exp[z/(f)x2] f p . /aN n 

' " ' """ -ig j (xp+px) exp(iTa) 


{Ai( 6 a; + a)} =- 


2tt 


dr exp ( i“^ ) exp 


X exp if ^ p^ exp {irbx). 


But, 


w I 


then, we have 

{Ai{bx + a)} = 


exp */ (^) exp (irbx) = exp i/ b 
exp (-if) exp [if (f) x"^] 


27r 


dr exp 1 ) exp 


exp (irbx ), 

exp(ira) 


(39) 

(40) 


Also 


exp 




X exp —ig ^ —^ {xp + px) exp (ibxr). 

(xp + px) exp (ibxr) = exp “5 (^) expjiexp -2^ (^) bxr'j , 
then the fractional Fourier transform of the Airy function can be written as 
5^„{Ai(5a;+ a)} = exp exp */(^) exp “5 (^) 

X— y dr exp y — ^ exp (ira) exp exp |i exp —2g(^ — 'j bxr'^ . 

Completing a cube binomial and changing the integration variable, we get 
S'a {Ai {bx + a)} = exp -i^ +if 


27r 


(41) 

(42) 


(43) 


(f) + *3^'^' (f) - (f) ""P [-25 (f )1 ^ - Kf) - (f) 

J drexp exp(ira)exp (^ir |exp -25 (f ) }) ' 


(44) 


Finally, recalling the integral representation of the Airy function (36) and the definitions of the functions / and g, 
we obtain 


{Ai( 6 a; + a)} =^/sec a exp < i 


Cul /^972 7*^ ^ 2 

— — + - tan a [ x — b a — xb sec a + — tan a 
2 2 \ 6 

54 

X b sec a + a —— tan a | . 


X Ai 


(45) 
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